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1. INTRODUCTION 


The primary objective of this research is to develop an efficient and robust trajec- 
tory optimization tool for the optimal ascent problem of the National AeroSpace Plane 
(NASP). Since the model of the aerospace plane is highly numerical data-driven, and very 
stringent flight path constraints must be imposed for the safety and operational reasons, 
it is felt that the existing trajectory optimization techniques either cannot handle this 
problem or do not offer a completely satisfactory solution in accuracy or efficiency. The 
issue of on-board guidance also motivates the research. Although the solution for the 
optimal trajectory may have to be generated off-line, it would be highly desired if some 
algorithm can be set up so that the actual flight information can be incorporated with 
the open-loop optimal solution in a feedback fashion to control the vehicle in real-time. 

This report is organized in the following order to summarize the completed work: 
Section 2 states the formualtion and models of the trajectory optimization problem. An 
inverse dynamics approach to the problem is introduced in Section 3. Optimal trajectories 
corresponding to various conditions and performance indeices are presented in Section 
4. Based on the accurate optimal solutions, a midcourse nonlinear feedback controller is 
develped in Section 5 to obtain significantly simplified yet very accurate optimal trajec- 
tories. The remarkabe performance of the inverse dynamics approach and the midcourse 
controller in guiding the aerospace plane in the presence of disturbances and uncertainties 
is demonstrated in Section 6. Section 7 discusses rocket-assisted ascent which maj< be 
beneficial when orbital altitude is high. The possibility of singular control of the rocket is 
examined. Feedback control laws for the rocket are derived by using the inverse dynam- 
ics concept. Finally, Section 8 recommends some future research aspects for the subject. 
Detailed discussions of some of the above topics can be found in Refs. [1-2]. 
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2. MODELS AND PROBLEM FORMULATION 


The model for the aerospace plane is adopted from Ref. [3], known as the Langley 
Accelerator”. It is a generic model of a hypersonic vehicle with winged-cone configuration. 
The present study is restricted to two-dimensional motion. The motion is controlled by 
propulsion throttle and angle of attack which in turn affects the aerodynamic forces. The 
propulsion system is assumed to comprise of only airbreathing engines for now. The 
thrust is modeled by 

T — Crq 


where q is the dynamic pressure, and Ct is the thrust coefficient. Ct a-nd the specific 
impulse I ap are given in tabular form as functions of Mach number, dynamic pressure 
and fuel-equivalence ratio, designated by <f>. Equivalence ratio of unity corresponds to 
maximum fuel efficiency, and values greater than unity give more thrust but use dispro- 
portionately more fuel. Assuming a spherical, nonrotating earth and Newtonian gravi- 
tational field with g = p/r 2 , the point mass equations of motion for the aerospace plane 
are 


dr 

dt 

dO 

dt 

dv 

dt 

d 7 

dt 

dm 

dt 


v sin 7 

v cos 7 
r 

T cos(a — e) — D 

m 

T sin(a — e) + L 

mv 


psmpf 


,v t 1 \ 

+ ( ) cos 7 

r vr 1 


T 

9o Isp 


( 1 ) 

( 2 ) 

(3) 

(4) 

(5) 


In above equations, r is the radius from the center of the earth to the vehicle, 9 the polar 
angle; v the velocity; 7 the flight path angle and m the total mass. T is the thrust and 
I S p is the specific impulse of the propulsion system, a denotes the angle of attack, e 
is the thrust vector angle and will be assumed to be constant or zero for airbreathing 
engines. L and D axe aerodynamic lift and drag, respectively. The atmospheric density 
is assumed to be an exponential function of altitude, 


P = Poe 


-P(r-ro) 
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where ro = 6378 km is the radius of the earth. The lift and drag coefficient Cl and Cd 
are also given as functions of Mach number and the angle of attack in tabulated data. 
For convenience of later application, Cl for the basic vehicle is accurately approximated 
by 

Cl = Clo<x + Clio: 3 ( 6 ) 

where Clo and Cli are functions of Mach number, obtained by least -squares fit to the 
tabulated data for given a, and then interpolated by cubic splines over Mach number. 
Figure 1 shows Clo and Cli as function of Mach number. Cl , C d and I S p are interpolated 
using a multi-dimensional table look-up scheme. 

The horizontal takeoff conditions from a runway at sea level are specified by initial 
velocity of Mach 0.5 and takeof mass of 133,809 kg (295,000 bf). Final orbit injection 
into a circular orbit at altitude h c is achevied by constraining the terminal altitude h(tf), 
velocity v(tf ) and flight path angle 7 (if) 


h(tf ) — he 
v(tf) = 


7 (</) = 0 


( 7 ) 


where the final time t j is free. Two operational constraints on the trajectory are 


9 — Qmax 

( 8 ) 

Q Qmax 

( 9 ) 


where Q is the convective heating rate at the stagnation point. The optimization is to 
be performed with respect to the control variable (f> and a so as to minimize the fuel 
consumption of the ascent trajectory, or, equivalently, 


max m{tf ) 

<t>, a 
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coefficients 



Figure 1 . Model of the lift: Ci — Ciqoi + Cl\oc 3 


4 



3. AN INVERSE DYNAMICS APPROACH 


A representative of the previous work is documented in Ref. [4] in which a singular 
perturbation method is utilized to analyze the near-optimal trajectory. Since no ana- 
lytical functional forms of the system parameters Co, Ct 3X1 d lap are readily available 
for our aerospace plane model, and the problem is strongly nonlinear, a direct method 
using nonlinear programming approach appears to be more advantageous for accurate 
solutions. But it was found that the optimization problem is very poorly conditioned. 
This is mainly because the aerospace plane flies at hypersonic speeds, and any slight 
variations of controls will cause large perturbations in the trajectory down stream. This 
high sensitivity problem makes even the conventional direct approach not convergent in 
which the control histories, a(t) and are directly parametrized. Early attempt to 

generate optimal trajectories succeeded only when no constraints (8) and (9) were en- 
forced, becuase the unconstrained trajectory accelerates through the dense atomsphere 
very rapidly (within 3 minutes, the altitude is already 42 km, and the velocity reaches 
Mach 20). Clearly the unconstrained trajectory will not be feasible. To circumvent this 
difficulty, an inverse dynamics approach is introduced. The inverse dynamic problem 
(IDP) can be formally stated as the following: 

Consider a dynamic system 


x = f (x(t),u(t),t) (10) 

with given initial conditions and terminal constraints 

x(t 0 ) = x 0 (11) 

S(x(t,),t/) = 0 (12) 

Find a control u (t) such that the solution of (10) with initial condition (11) satisfies (12) 
and the given algebraic constraint 

g(x(<),c(<),<) = 0, <e[*o,</] (I 3 ) 

where g : R” x R m x R — > R* is sufficiently differentiable, c (t) E R m for to < t < tf is 
a given smooth function. c(t) usually represents the desired output and (13) specifies the 
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output relationship. By repeatedly differentiating each component of (9) till u appears 
explicitly, we have additional constraints 

G(x(f), u (<), c (<), c (t ), ..., t) = 0 (14) 


Equations (13) and (14) are the constraints on the state variables and controls. The 
existence of solution to such an inverse dynamic problem can be guaranteed under certain 
conditions [5]. 

Most work in this area has so far centered on finding the required control for a given 
c (t). We have extended this idea to trajectory optimization. In a standard IDP, the 
original nonlinear system needs to be transformed into a system linearly dependent on 
the control vector [6] . This transformation facilitates a general methodology to solve for 
u. But for the current problem, our following treatment is more advantageous. To apply 
the idea of inverse dynamics more effectively, we first change the independent variable 
from t to 6. The system equations (l)-(5) with e = 0 now become 


dr 

— = r tan 7 

dt _ r 
dS v cos 7 
dv ,T cos a — D 
d9 m 


^sin7 


v cos 7 


dj ,Tsma + L , ,v n ^ x r 

— = ( V ( -)cos7) 

d6 mv r vr 2 v cos 7 


mv 

dm _ T r 
d9 golsp v cos 7 


(15) 

(16) 

( 17 ) 

(18) 
(19) 


Analogous to (13), we define 


g = r(0) - e(«) = 0 


( 20 ) 


In (20) c(B ) is a sufficiently smooth function which represents a specified altitude history. 
Differentiating (22) once with respect to 6 gives 



( 21 ) 


Differentiating (21) once again, 


, N .... 9 ,ucos' J 7 ,v n . Tsino. 

L(a ) = mu (c - r tan 2 7) ■= ( ^)cos7 ] 

K ' lv r 2 r vr 2 mv 


( 22 ) 
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The prime in (21) and (22) stands for differentiation with respect to 0. Equations (21) 
determines the required 7 for r to follow c(0). Equation (22) provides the necessary lift 
control. For a specified thrust level T and current values of the state variables, Equation 
(22) constitutes an algebraic equation for a. With a solved from (22), Equations (17) 
and (19) can be integrated for v and m at next instant. So the solution of the system as 
well as the value of the performance index (fuel consumption, for instance) is completely 
determined by the choice of the pair of command altitude c and fuel-equivalence ratio 
If we choose to represent c($) and <^($) by certain smooth parametrized functions, 
the optimization problem reduces to a parameter optimization problem in which the best 
c and <{> histories are iteratively sought through solving a sequence of inverse dynamic 
problems. Note that in the parametrization of c, one can always choose the boundary 
conditions 

c(0 f ) = r o + h e (23) 

c\0 f ) = 0 (24) 

Then the first two of the three terminal constraints in (7) are automatically satisfied 
according to (20) and (21), leaving only the constraint on Vf to be met. With Cl 
represented by (6), a can be solved from (22) very effectively by Newton iterations with 
an accuracy of 10 -6 frequently after only one iteration. Thus the computation of a does 
not pose extra burden since r and 7 equations do not need to be integrated as a result. 

With this inverse dynamics approach the trajectory is under more direct control of 
the parametrization process. Consequently, the sensitivity of the optimization problem is 
greatly reduced. In fact, with minimum efforts, one can easily construct various feasible 
trajectories that satisfy the terminal conditions (7) and state constraints (8) and (9) by 
choosing c(0 ) and <f>{9). This feature is not only essential to the success of the trajectory 
optimization, but may also be useful for quick design of hypersonic cruising trajectories 
(not necessarily optimal in any sense). Compared to the collocation method in Ref. [7], 
the current approach retains the merits such as better conditioning of the problem and 
robustness of the algorithm, while it cuts down the dimension of parametrization consid- 
erably for an accurate solution. Fewer optimization parameters directly contribute to a 
faster convergence. In addition, the inverse dynamics approach permits guidance com- 
mands to respond adaptively to perturbations, which will be discussed in later sections. 
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4. OPTIMAL SOLUTIONS 


In order to gather sufficient information on the optimal trajectories, extensive runs 
have been performed to generate fuel-optimal trajectories. First, only constraint (8) is 
considered. The constraint (8) is handled by introducing a new state w q 


W q (0) =0, w q = 


f — (? — Qmax)i 

u 


9 > ?maz! 

9 — Qmax • 


(25) 


Constraint (8) is then equivalent to the terminal constraint 


W q (9f) > 0 


(26) 


A similar transform can be done for constraint (9). A squential quadratic programming 
algorithm [8] was used to solve the resulting constrained nonlinear programming problem. 
A parametric study on q m ax was conducted. Figure 2 depicts variations of q for different 
values of q max ■ In Fig. 3 are fuel-equivalence ratios for three trajectories. Figures 4 and 5 
contain typical flight path angle and angle of attack histories and ascent altitude history. 
The wiggles on the g-histories in Fig. 2 are more likely a result of the finite- dimension 
approximation of the original optimal control problem than anything physically signifi- 
cant. The final masses for the optimal trajectories corresponding to different constraint 
levels axe listed in Table 1. 


Table 1. Summary of Optimal Solutions 

q — qmax 


Qmax (psf) 

OO 

5000 

4000 

3000 

2000 

™(</) ( k s) 

69859 

68781 

68525 

68101 

66351 


General features of the optimal trajectories include: 

(1) The trajectory typically divides into three stages: a quick climb out following 
takeoff, featuring large flight path angle and angle of attack; a relatively long mid-course 
cruise with small flight path angle and angle of attack, during which the constraint (8) 
is active; a final zoom into orbit. 

(2) The optimal fuel-equivalence ratio stays almost unity during the midcourse cruise. 
Some modulation may be needed in the initial climb out to avoid violation of the con- 
straint (8). 
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(3) When g max > 2500 psf, the penalty of the constraint (8) on the fuel consumption 
is not very significant. The optimal trajectory simply takes longer flight time along the 
boundary q = q m ax to accelerate to the required orbital velocity. When q m ax is below 
2500 psf, the drag loss during the long cruise becomes significant, resulting in considerable 
increase in fuel consumption. 

The observation (3) above prompts the question of what is the smallest q max below 
which the NASP cannot achieve the orbital velocity with all the fuel it carries. The 
performance index this time should be 

J = max q(t) 

to <t<tj 

The corresponding optimal control problem is the so called minimax problem. The fol- 
lowing transform converts this nonclassical optimal control problem into a parameter 
optimization problem, hence the present approach is still applicable. 

min Qmax 


Q £ Qmax 
7Yl(t f ) ^ rri m i n 

subject to system equations (l)-(5) and boundary conditions (7), where ra mm is the mass 
of the aerospace plane excluding fuel. For our model, m mtn = 63,503 kg (140,000 bf). 
qmax now is treated as an optimization parameter. For injection into the same orbit, 
the minimum q ma x turns out to be 1499 psf with all fuel consumed (70,306 kg). The 
variation of q in this case is also plotted in Fig. 2. 
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(jsd) b 



Figure 2. Variations of dynamic pressure q 
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fuel-eqv. ratio 



Figure 3. Variations of fuel-equi valence ratio <f> 
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Figure 4. Flight path angle 7 and angle of attack a along the optimal 
trajectory q < 2000 psf. 
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Figure 5. Optimal ascent altitude history ( q < 2000 psf). 
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5. MIDCOURSE FEEDBACK CONTROLLER 
AND SIMPLIFIED SOLUTIONS 

As we have observed, the prominent features of the midcourse include slowly- varying 
'y, small a, fuel equivalence ratio of nearly unity, and observance of q — qmax • On the 
basis of these observations, we can develop a nonlinear feedback midcourse controller that 
maintains these features of the trajectory. 

We set (j) = 1 after initial climb out. Since 7 & 0 and a is small in the midcourse, 
letting 7 = 0 and neglecting T sin a in Eq. (4) lead to 

L = m(-^- - — )cos 7 (27) 

H r 

(27) is essentially the same as the reduced solution in Re.f [4] in which it is derived based 
on a singular perturbation method. In light of our optimal solution, the validity of the 
time scale decomposition is confirmed for the midcourse. Nonetheless, (27) is not a valid 
approximation for the initial climb during which both 7 and 7 have large values. 

Controller (27) does not have the mechanism to enforce state inequality constraints. 
However, we have observed that the optimal trajectory climbs on boundary q = q m ax an d 
Q = Q m ax in sequence during the midcourse. To renforce the constraints on the basis of 

(27) , a midcourse feedback controller 

_ j m{n/r 2 - + h(q - q max ) + hq, if Q < Qmax 

G= -v^M^^ + hiQ-Q^ + hQ, if Q> Qmax 

is proposed. The feedback gains k t ’s are properly chosen constants. In fact, controller 

(28) can be used to track any hypersurface F — F re f = 0 when F is a given function of 
the altitude and some other state variables. If F is a monotonic function of the altitude, 
which is the case for q and Q, constant coefficients ki's will suffice. 

With the aid of controller (28), the simplified trajectory is planned in the following 
sequence: For the initial climb 9 € [0, $i), both 4>{ff) and c(0) are parametrized by cubic 
splines as before. For the midcourse 9 € [# 1 , $ 2 ), let <f> = 1, and a be given by (28). The 
final zoom part 9 6 [^2 >$/] is characterized by letting <j> = 1 and 

c(9) = a + b9 + d9 2 + e9 3 (29) 

Equation (29) determines the angle of attack control for the final zoom via the inverse 
dynamics. 9\, 92, 9f , c($) and <^(^), 9 E [0,$i], are to be optimized. The coefficients in 


14 


(29) are determined by the continuity conditions c{9 2 ) — r(# 2 )j C '(Q 2 ) — r(^ 2 )i ;an ^ 2 ? 
the terminal conditions (23) and (24). 

Planned in such way, the dimension of the optimization parameter vector is reduced 
almost by two thirds, down from over 30 to 12, which contributes to a much faster con- 
vergence of the optimization process. More significantly, the conditioning of the problem 
is further improved to such an extent that the rather difficult original optimization prob- 
lem now becomes an easy routine operation. The fuel consumption is further reduced 
because the midcourse controller tracks the boundary q = q m ax better than the control 
a obtained through finite-dimension parametrization, hence the efficiency of the the air- 
breathing engine is better. For instance, the constrained solution with q < 95, 760 N/m 2 
yields a final mass of 67,112 kg, greater than the previously obtained value of 66,351 kg in 
Table 1. Figure 6 shows the ^-history and Q- history along the optimal trajectory under 
the controller (28) with q < 95,760 N/m 2 and Q < 800 Watt/cm 2 , where the heating 
rate Q is modeled by 

Q = (4.919 x 10“ 8 )/ V o (30) 

In (30) density p is in kg/m 3 and velocity v is in m/sec, and (30) corresponds to equilib- 
rium conditions on the surface of a wing leading edge 10 cm in radius [4]. Shown in Fig. 
7 is the ascent altitude history for the same trajectory. 

In closing this section, we conclude that control law (28) not only greatly improves 
the process of off-line generation of the optimal solution, but is also particularly attractive 
for real-time onboard applications, because the controller demands no intensive calcula- 
tions for implementation. The simplified optimal solution provides an fairly efficient and 
reliable way of investigating the optimal trajectories of a very difficult problem. Because 
of the feedback nature of the midcourse controller, even under the off-nominal flight 
conditions the controller is still able to prevent sizable violation of constraints (8) and 

(9). 
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q(psf). Q(Watt/cm**2) 
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Figure 7. Optimal ascent altitude history ( q < 2000 psf, and Q < 800 
W/cm 2 ). 
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6. ASCENT GUIDANCE 


In addition to the ever exsiting environmental disturbances, the unprecedented na- 
ture of the NASP increases the likelihood of system uncertaintities and modeling inaccu- 
racy. All these factors make it imperative that the NASP have a robust guidance system 
that can direct the vehicle to fulfil its mission in the presence of the disturbances and 
perturbations. The inverse dynamics approach introduced in Section 3 renders a powerful 
tool in this regard. 

Once a nominal optimal trajectory for a flight vehicle is designed, there are two types 
of popular guidance strategies. In the first strategy corrective actions based on devia- 
tions of the actual trajectory from the nominal are taken to restore the flight path. The 
second strategy recalculates a new trajectory on-line based on current state if deviations 
occur, using techniques such as neighbouring extremal control and singular perturba- 
tions. As we shall see, the inverse dynamics approach in conjunction with appropriate 
feedback compensation offers a blend of the above two strategies, and appears to offer 
both robustness and suitability for real-time environment 

Suppose that the equivalence ratio <f>* and command altitude c have been obtained 
off-line for an optimal ascent trajectory of the NASP. The generation of the corresponding 
angle of attack command a from 


L(a) — mv[(c*" — r tarn 7 ) 


2 rcos J 7 ,v r . 

,2 -A L _ ( -) cos 7 - 


±L 

vr 2 


T sin a , 


mv 


(31) 


only requires minimum computation and should be easily carried out onboard. The 
values of the state variables in (31) are the actual values instead of nominal ones, a so 
generated is adaptive to perturbations and disturbances when they are known. The idea 
is equivalent to solving the inverse dynamic problem for a different system to generate 
the same output. For instance, if the atmospheric density has fluctuations, as it always 
will happen, and they are measured instantly, the terms involving p in Eq. (31) can 
be assigned the correct value, a so computed would result in dy j d9 = dy*/d9, then 
•y( 0 ) = y*(9) and r(9) = r*(9) if there are no initial errors, where the quantities with 
asterisk are the nominal values. Nonetheless, perturbations and disturbances are not 
always known to the onboard computer. Additional feedback terms may be added to 
compensate the inaccuracy in more general situations. A such candidate is 


a = ci — kft(h - h *) 


(32) 
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where fc/, > 0 is guidance constant. When perturbations and/or disturbances are known, 
a reduces to d since h = h*. The deviations in velocity are very effectively compensated 
by fuel-equivalence ratio 

<f> = (f>* — k v (v — v*), k v > 0 (33) 

To test the effectiveness of (32) and (33), the simplified fuel-optimal trajectory with 
q < 95,760 N/m 2 is chosen to be the nominal. Equations (32) and (33) are applied 
in the initial climb and final zoom. The midcourse is still controlled by (28) without 
modification, which assures that no excessive violation of the state constraints would 
occur along the perturbed trajectories. Some typical sources of errors are examined in 
the following. 

Density Fluctuations 

Assume that the actual atmospheric density is varying according to 

p = (1 + 0.2 sin ^^)p* (34) 

In (34) p* is the nominal density, h is altitude in km. The maximum fluctuation of 20% 
as given by (34) is a modest variation of the density compared with the flight data of the 
Space Shuttle [9]. 

(a) The variation is measured and known to the onboard computer. By the above 
arguements, h and 7 will follow their nominal histories exactly. With (33) the final error 
in velocity was only A vj = 0.97 m/sec. 

(b) Suppose that the density variation is not known to the onboard computer. Figure 
8 exhibits the comparison of the nominal trajectory with two perturbed trajectories: (1) 
Perturbed trajectory I used the nominal density to calculate a and the controls were 
obtained from (32) and (33). Final orbital insertion errors were Ay/ = -1.66 m/sec, 
A h f = 0.747 km and A7 / = 0.37°. (2) Perturbed trajectory II was flown with the 
preprogrammed nominal optimal a* and (f >* . The open-loop controls resulted in a crash. 

Lift Coefficient Variation 

Suppose that a uniform inaccuracy in Cl exists such that Cl only has 90% of its 
nominal value 

C L = 0.9Cl (35) 
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(a) If the inaccuracy can be identified in flight, ot can be calculated using the true 
values of Cl- The only terminal error was in v which was reduced to A Vf = —0.72 m/sec 
by (33). 

(b) More likely the true values of Cl will not be exactly known, a will be computed 
with the nominal value C* L . Two perturbed trajectories, again denoted as perturbed I and 
II as in above, are plotted in Fig. 9. Trajectory perturbed I which was controlled by (32) 
and (33) nearly followed the nominal. The difference is virtually disernible to the scale 
of the plot. Final errors were A Vf = —0.99 m/sec, A hf = 0.17 km and A 7 / = 0.021°. 
The open-loop nominal controls failed to steer perturbed II to achieve the final, orbit. In 
fact, the two dives along trajectory perturbed II generated huge peak dynamic pressure 
(>2,000,000 N/m 2 ) that would have already crushed the vehicle. 

Thrust Inaccuracy 

Another source of perturbations is the model inaccuracy of the airbreathing propul- 
sion system. Suppose that a 5% loss in actual thrust level attributed to the model 
inaccuracy exists: 

T = 0.95 T* (36) 

The largest error in this case will be in velocity, since the flight path is mostly con- 
trolled by the angle of attack. The simple compensation (40) helps reduce the velocity 
error remarkably, resulting in A Vf = —0.73 m/sec while A vj = 158.3 m/sec without 

compensation in <f>. 
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Figure 8. Comparison of trajectories in the presence of atmospheric den- 
sity fluctuation: nominal=optimal trajectory ( q < 2000 psf); perturbed 
I=trajectory flown with guidance laws (27), (28) and nudcourse controller 
(26); perturbed II=trajectory flown with preprogrammed nominal controls 
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Figure 9. Comparison of trajectories in the presence of lift modeling inac- 
curacy: nominal=optimal trajectory (q < 2000 psf); perturbed I=trajectory 
flown with guidance laws (27), (28) and midcourse controller (26); perturbed 
II=trajectory flown with preprogrammed nominal controls 
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7. ROCKET-ASSISTED ASCENT 


7.1 Singular Control Analysis 


When the orbital altitude is higher, during the final push the airbreathing propulsion 
system, supersonic ramjets (SCRAM jets), may become ineffective as the atmosphere gets 
thinner. A rocket will be needed. Calise et al show that the SCRAM jets may still stay 
on for the remaining flight even if they are not effective [10]. Practically, the SCRAM 
jets may be cut off at an optimal point. Assuming a throttable rocket, the remaining 
rocket-assisted trajectory could consist of a combination of coasting, singular thrust and 
full-throttle arcs. We shall examine the possibility of singular arc first. Let a = 0 after 
the SCRAM jets are turned off. By the standard optimal control theory [11] and system 
(l)-(5), we have the Hamiltonian 


T cose ,D . sin 7 N Tsine , _ ,v n N „ 

H = p r V Sin 7 + Pw Pv(— + rr „, +P7( r “ yj.2 ) cos 7 Pm 


771 


771 


mv 


90 1 


sp 


= 0 

(37) 


The switching function is 

. mvpm 

S = vp v COS £ — Py Sin £ Z 

ffo *sp 

where for the rocket, I sp is considered constant. The optimal thrust is given by 

r 0, S < 0; 

T* = \ 0 < T* < T m ax, S = 0; 
l Tmax, S> 0. 


(38) 


(39) 


The first case in (39) corresponds to a coasting arc, second to a singular arc. The optimal 
thrust angle is derived from dH j de = 0: 


tane* 


P 7 
VPv 


(40) 


By the minimum principle [10], the matrix 


( d 2 H/de 2 d 2 H / dedT \ 
\d 2 H/dTde d 2 H/dT 2 ) 


(41) 


must be positive semidefinite along an optimal solution. It is obvious by (37) and (40) 
that d 2 H/dT 2 = 0 and d 2 H/dedT = d 2 H/dTde = 0. On singular arc one can show by 
using 5 = 0 that 

d 2 H _ mp m v 
de 2 go hp 
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Since p m represents the sensitivity of the minimization performance index J = -m(tf) 
with respect to variation of m, p m should always be negative. In fact, if drag D = 0, 
one can arrive at, with the aid of adjoint equations p = —dH/dx and transversality 
conditions, the relationship 

p m (t) = -ZM (43) 

m(t) 

In conclusion, we see that d 2 H/de 2 > 0, or equivalently, H uu is positive semi-definite. 
The necessary conditions for optimality do not exclude the possibility of singular arc. 
Other tests for the optimality of singular control (e.g. Ref. [12]) are difficult to apply 
here due to the complexity of the system. But given the observation that at higher 
altitudes where aerodynamic drag is small, a sustained thrusting arc can be beneficial 
to the buildup of the insertion velocity, singular arc may exist. However, analysis in the 
Appendix asserts that if the final orbit is circular, the optimal arc immediately before 
orbital insertion must be a full-throttle one. 

On singular arc, it is shown in the Appendix that z = tan(e/2) satisfies 

Az 4 +Bz 3 +Cz 2 + Ez + F = 0 (44) 


where A, B , C, E and F are functions of the state variables. In particular, if |e| << 1, 

we have , „ _ . r . , 

sin 7 [Dr 2 (v + gplsp) — sm 7j 


£ ~ _F 

tan g ~ £ 2(mpgoIsp sin 2-y + r^vD cos 7) 


(45) 


The throttle on a singular arc can in principle be obtained in a similar way, only the 
expression turns out to be excessively complicated. 

Within the framework of nonlinear programming approach, the possibility of coasting 
and singular thrust arcs is conveniently handled. The complete ascent trajectory with 
rocket assistance is parametrized in the following sequence: airbreathing-engine powered 
ascent portion as before; a coasting arc following the cutoff of the SCRAM jets with 
a = 0; rocket- assisted flight to orbital insertion. The optimization parameters include all 
those used in previous sections plus time durations of coasting and rocket-assisted flight, 
as well as those that specify the rocket throttle and e programs. 


7.2 Feedback Control Laws via Inverse Dynamics 

Despite that the analysis suggests possible combination of singular/full-throttle arcs, 
it was found that the fuel-consumption is not sensitive to the rocket throttle program at 
all. The reason is that after a long coast the rocket is turned on at approximately the 


24 


I 


orbital altitude. The optimal |e| then is almost zero and the major role of the rocket 
is to push the velocity to orbital speed. Since drag D is very small at that altitude, 
the situation is similar to a rocket with fixed £ in vacuum— the velocity increament is 
only dependent on the rocket fuel expenditure, independent of the throttle program. 
Numerical results show that allowing throttle to vary between [0,T moi ] yields negligible 
improvement in the fuel-consumption (only few kilograms of fuel saving) as compared 
with the case of full-throttle. Nevertheless, using intermediate thrust level of a throttable 
rocket can be rather advantageous in achieving orbital insertion accurately. To show this, 
we again adopt the inverse dynamics approach. Let the command altitude c(9) and a 
desired history of velocity u($) be specified for the rocket-assisted portion of the trajectory. 
In a similar way as Eq. (22) is derived, we have 


. ,v n \ 

T sme = 5- mu cos 7 - 

r vr z 


mv 2 cos 3 7 „ 


2 {c" — r tan 2 7) = U 


T cos £ = 


mu cos 7 , usin7 
~v — 7, ” 


+ D = V 


(46) 

(47) 


Two resulting feedback control laws axe 


v (48) 

max 

tan £ = ^ (49) 

where 77 represents the throttle of the rocket, and T max the maximum available thrust. 
Instead of optimizing c(9) and u($), we find it sufficient to let 


c(9) = a + b9 + d9 2 + etf 3 (50) 

v(9) = n + p9 (51) 

With the rocket thrusting duration being an optimization parameter, the coefficients a, 
&, d and e in (50) are determined by the continuities of r and 7 at the instant 9 r when 
the rocket is turned on, and by Eqs. (23) and (24) for the first two terminal constraints 
in Eqs. (7). Coefficients n and p are defined by the continuity of u at 9 r and the final 
constraint in (7). Since these calculations can be easily done onboard, the near-optimal 
parametrizations (50) and (51) have a distinct advantage: the orbital insertion conditions 
(7) will remain satisfied even if the actual trajectory deviates from the nominal at 9 r . 
This is because (50) and (51) always lead the trajectory from the states at 9 r , whatever 
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the values are, to the target point defined by Eqs. (7) when the coefficients are calculated 
with the actual values of the states at 9 r . The corresponding inverse controls (48) and 
(49) are hence disturbance-accommodating, provided that they are not in violation of the 
control constraints |e| < £max and tj < 1. 

A typical near-optimal rocket-assisted trajectory for circular orbital insertion at hf — 
92.6 km (50 nm) subject to constraints (8) and (9) is plotted in Fig. 10. The points of 
interest are marked. In particular, the SCRAM jets are found to be turned off at 70 km. 
Following a coast of 452 seconds, a throttable rocket with T max = 266, 893 N (60,000 lb) 
and I sp = 440 seconds starts firing at 92.3 km with an intermediate throttle of about 
0.21. Orbital insertion is achevied 88.5 seconds after the rocket ignition. The final mass 
is 64,931 kg, which is comparable to m/ = 65, 370 previously obtained for hf = 55 km 
without rocket- assistance. It should be noted that the coasting arc saves fuel significantly 
(more than 3,000 kg) by allowing at no cost the change of kinetic energy generated by 
the airbreathing propulsion to potential energy at higher altitude. 

7.3 Guaranteed Orbital Insertion 

The guarantee of accurate orbital insertion with the inverse control of the rocket in 
the presence of disturbaces and perturbations is demonstrated by introducing atmospheric 
density fluctuations. Vertical variation (34) and the following longitudinal variation 

p= (1 + 0.2 sin ^q)p* ( 52 ) 

are used, where in (52) d = r 0 B (km) is the down range. Assuming that the actual p is not 
known to the onboard computer, the airbreathing propulsion portion of the trajectory is 
controlled by the midcourse controller (28) and the guidance laws (32) and (33). These 
measures have already greatly reduced the deviations of the perturbed trajectories, which 
prevents the trajectories from entering unrecoverable regions of the rocket. Then after the 
coast the inverse rocket control laws (48) and (49) adaptively steer the aerospace plane 
to an accurate orbital insertion. Figure 11 contains the nominal and the two perturbed 
trajectories. The final masses of the perturbed trajectories are 64,635 kg and 64,202 kg, 
respectively. The longitudinal p fluctuation is more harmful to trajectory control than 
the vertical deviation. This is no supprise if one compares the down range the aerospace 
plane travels (about 14,300 km) with the altitude gained (92.6 km). The oscillations of 
altitude along the perturbed trajectory II in Fig. 11 suggest that the midcourse controller 
(28) is at work, tracking the q and Q-constraints accurately even in the presence of the 
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sinusoidal perturbation of the atmospheric density. The variation of the corresponding a 
is shown in Fig. 12. Figure 13 depicts the histories of q and Q in this perturbed situation. 
The nominal and perturbed rocket throttle settings r) and thrust angles e are shown in 
Figs. 14 and 15. The effect of controls adapting to disturbances are clearly seen. 
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Figure 10. Typical rocket-assisted optimal ascent 
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Figure 11. Comparison of rocket-assisted trajectories: nominal =optimal trajectory in Fig. 
10; perturbed 1= trajectory subject to vertical atomspheric density fluctuations: perturbed 
II=trajectory subject to horizontal atomspheric density fluctuations. 
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alpha (deg) 



Figure 12. Variation of angle of attack along the trajectory (perturbed II) 
subject to horizontal sinusoidal atmospheric density fluctuations. 
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rocket throttle 



Figure 14 . Comparison of rocket throttle settings along rocket-assisted 
trajectories (perturbed I and II are defined the same as in Fig. 11). 



thrust angle (deg) 



dimensionless time 


Figure 15. Comparison of rocket thrust angles along rocket-assisted trajec 
tories (perturbed I and II are defined the same as in Fig. 11). 
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8. CONCLUSIONS AND FUTURE RESEARCH 


The trajectory optimization problem for the National AeroSpace Plane ha s been for- 
mulated as an inverse dynamic optimization problem. Extensive numerical expreiments 
have been conducted. The approach has proven effective in solving this otherwise very 
difficult problem with relative ease. Generalization of similar approach to some other 
difficult optimization problems may be worthwhile. Both classical performance index - 
minimum-fuel and nonclassical performance index - minimax dynamic pressure are con- 
sidered. The results lead to a better understanding of the characteristics of hypersonic 
flight. A nonlinear feedback midcourse controller suitable for onboard implementation 
is proposed, which also significantly further simplifies and improves the solution. Ro- 
bust ascent guidance is obtained by using combination of feedback compensation and 
onboard generation of control through the inverse dynamics approach. Optimal rocket- 
assisted trajectories axe investigated. The pattern of the optimal tajectory is found to 
be airbreather-powered ascent + coast + rocket-assist. The last arc immdediately before 
insertion onto a circular orbit must be a full-throttle arc. Inverse rocket control laws 
are developed that together with the ascent guidance scheme guarantee accurate orbital 
insertion even in the presence of disturbances and system uncertaintities. This work has 
also prompted the following thoughts that warrant further investigation: 

1. Analysis of Constrained Trajectories 

As found in above sections, the dominant portion of the optimal ascent trajectcory 
lies on the boundaries of the dynamic pressure and heating rate constraint. It is believed 
that accurate approximate analytical solution of time to the constrained portion of the 
trajectory is possible. An enlightening treatment for a special case is presented in Ref. 
[13] in which the constrained dynamic system is shown to be a two-time-scale system and 
asymptotic analytical solution is obtained. The idea begins with a careful examination 
of the constrained flight of an aerospace vehicle 

x(t) = f(x,u) (53) 

C(h,p(h),v) = 0 (54) 

where Eq. (53) is the state equation and x is the state vector which typically consists of 
h (altitude), 8 (polar angle), 7 (flight pathe angle), v (speed) and m (mass), u is the con- 
trol. Equation (54) is a nonlinear algebraic flight path constraint, p in (2) represents the 
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atmospheric density. By properly choosing a set of nondimensional variables and study- 
ing the predominant factors of the constrained motion, it may be possible to reveal that 
a natural two-time-scale property exists so that the flight path angle dynamics can be 
considered “fast” as compared to the altitude dynamics. With the aid of singular pertur- 
bation theory, it may be possible to obtain approximate analytical asymptotic solutions 
of time. Then the trajectory optimization and guidance problem will be significantly 
simplified. 

2. Application of Simulated Annealing Algorithm 

While the inverse dynamics approach has been quite successful, it requires the prob- 
lem to be sufficiently smooth. The smoothness is achieved by using analytical fitting to 
the key aerodynamic coefficients. This may or may not be done in a reasonably efficient 
way, depending on the problem. In our previous endeavors this was fortunately the case. 
If a complete model of the aerospace plane is totally based on linear interpolations of 
tabulated data, which is most likely for a realistic vehicle, the trajectory optimization 
problem is inherently nonsmooth, although continuous. The nonsmoothness can cause 
serious convergence problems if optimization algorithms based on gradients are used. A 
genetic algorithm, a nongradient type method, has been tested on a similar problem . We 
will propose to use another type of nongradient method, continuous simulated annealing 
aigoritm [14], on the problem. The simulated annealing algorithm is a stochastic method 
that promises to find the global optimum in probability. In addition to strong robust- 
ness, it has other merits such as easy implementation and being suitable for large scale 
problems in that increasing the number of optimization parameters results in very small 
increase in computational time, which seems to be what genetic algorithms and other 
nongradient methods lack. The simulated annealing algorithm has not been applied in 
trajectory optimization so far. The application, if successful, will not only provide a badly 
needed alternative to the trajectory optimization problem for the aerospace plane, but 
will also establish a promising approach for other optimzation problems of nonsmooth 
dynamic systems. 

3. 3-D Maneuvers 

Our preceding work has revealed that the two constraints (8) and (9) severly restrict 
the flight envelop of the NASP in 2-D flight. The fuel-optimal ascent trajectory does 
not enjoy sinificant fuel saving as compared to a feasible trajectory that satisfies all the 
constraints. 3-D maneuvers may relax some of the restriction, thus offer better fuel 
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economy. Admittedly, 3-D optimal solutions will be even more difficult than the 2-D 
case. The investigation toward this direction necessitates a more robust optimization 
algorithm. The simulated annealing algorithm may help in this regard. 

4. Adaptive Guidance 

It has been observed that a set values for the guidance parameters that work very 
well for certain disturbances may not be the best choice for other disturbances. In fact, 
the ranges for the best choices of the guidance parameters for different disturbances that 
may be encountered in the flight of the aerospace plane may have no overlap. This 
suggests that an adaptive gain update scheme may remarkably improve the performance 
of the guidance system. Since the system dynamics is highly nonlinear, some new adaptive 
algorithm needs to be developed. The application of neural network to form a closed-loop 
adaptive mechanism for determining guidance parameters may prove attractive. 
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APPENDIX 


Singular Rocket Thrust Control 


The adjoint system to Eqs. (1), (3)-(5) is 


P„dD 2^si „ 7 . 2, 

m dr r 3 r 2 vr 3 

p m dD p v T sine , p , 

p v — — p r sin 7 2 Pyi. ~ d 2 2 ) cos ^ 

yv y m dv mv 2 r v l r l 

/icos 7 

P 7 — -p r U COS 7 + r 2 

(A 1) 
(A 2) 

p \ • 

+ P 7 ( t ) sin 7 

r t>r z 

p v T cos £ p 7 T sine 
— m 2 m 2 t> m 2 

(A3) 
(A 4) 

On a singular arc, the switching function 5 = 0 


mv Pm n 

S — vp v cos e — sin s = U 

M5) 

dH/de = 0 leads to 

i>p v sin e + p 7 cos £ = 0 

(46) 

The semi-positiveness of H uu in Eq. (46) depends on the sign of 


dH 2 

— — = - V p v cos e + p 7 sin e 
de l 

(47) 

With Eqs. (A 5 ) and (A 6 ), Eq. (A7) is equivalent to 


dH 2 mp m v 

de 2 ~~ gohp 

(48) 

which is Eq. (42). As pointed in the paper, p m < 0. In particular, if D = 
arc is singular, using Eqs. (A4)-(A6) and p m (tf) = —1 results in 

0 and the final 

m(t)p m (t) = constant = —m(if) 


For minimum-fuel problem with free end time, the Hamiltonian H - 

= 0 (Eq. (37)). 

On a singular arc, we also have 


. D sin 7 . 
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r vr l 

(49) 
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Differentiating (A9) twice and after much algebraic operation, we obtain for z = tan(e/2) 


Az* + Bz 3 +Cz 2 + Ez + F = 0 


(A10) 


where 
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When the final orbit is circular, the above conditions can be used to establish that 
singular arc is not fuel- optimum solution for the portion of the trajectory immediately 
before orbital insertion. Wc show this by contradiction. Assuming that the last optimal 
arc is singular, substituting orbital insertion conditions (7) in Eq. (A9) gives 


Pv(tf} D _ Q f ^ = Q s i nce d ^ 0 

m{tf) 


(A 12 ) 


Also combining conditions (9), (A10) and (All) produces 


e(tf) = 0 (^413) 

Then Eqs. (A5), (A 12 ) and A(13) give rise to 

Pm{tf ) = 0 ("414) 

which is contradictory to the transversality condition p m (f/) — — 1 - So singular arc is 
not optimal. On the other hand, since at the orbital altitude drag D is very small, a 
coast arc cannot achieve circular orbital insertion for any reasonable length of time (In 
particular, if D = 0, circular orbital insertion is impossible by coasting). In conclusion, 
the final arc of a fuel-optimal trajectory must be full-throttle. 
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